3.1 Fitting a model: explicit calculation

The maximum likelihood equations for generalized linear models are non-linear and do not have closed form solutions, hence we will develop a numerical approach. We’ll do this here from first principles.

So that we have a clear understanding of what is going on, we will focus on the simplest example: an intercept-only model (so simple in fact that we can write down the exact solution). We’ll then extend the ideas to our original two-parameter model.

\[E(Y_i) = \lambda. \]

With data \(y_1,y_2,\ldots y_n\) our likelihood is

\[L(\lambda)= \prod_{i=1}^n \exp(-\lambda)\frac{\lambda^{y_i}}{y_i!} =\exp(-n\lambda)\frac{\lambda^{\sum_{i=1}^ny_i}}{\prod_{i=1}^n y_i!} \]

then the log likelihood is

\[l(\lambda) = -n\lambda + \log(\lambda) \sum_{i=1}^n y_i -\sum_{i=1}^n \log(y_i!). \]

We can easily maximize this function by looking at the derivative

\[ \frac{dl}{d\lambda} = -n + \frac{\sum_{i=1}^n y_i}{\lambda} \]

so if the derivative is zero, must have \[\hat{\lambda} = \frac{\sum_{i=1}^n y_i}{n}. \]

The second derivative is

\[ \frac{d^2l}{d\lambda^2} = -\frac{\sum_{i=1}^n y_i}{\lambda^2} <0, \]

so we clearly have a maximum for

\[ \frac{d^2l}{d\lambda^2}\vert_{\lambda = \hat{\lambda}} = -\frac{n}{\bar{y}} <0. \]

But suppose that we could not explicitly solve for \(\widehat{\lambda}\). We could instead proceed numerically, e.g. by using Newton’s method.

We start by giving the log likelihood function for the intercept-only model

llk_pois <- function(x) sum(dpois(y,
                                  lambda = x,
                                  log = TRUE ))

Let’s simulate some data from the intercept-only model and plot the log likelihood as a function of its parameter.

lambda <- 3
y <- rpois(n, lambda = lambda)
lambda_seq <- seq(.1,6, length.out = 1000)

llk_curve <- sapply(lambda_seq,
             FUN = llk_pois )

plot(lambda_seq,
     llk_curve,
     type= "l",
     xlab = expression(lambda),
     ylab="log lk")

We will code up a method for finding the maximum point of this curve numerically.

We do this by successively maximizing approximations to the curve.

3.2 Fitting the model with Newton’s method

We’ll see how Newton’s method converges rapidly to the maximum.

We expand the log likelihood as a Taylor series to second order. For small \(h\),

\[ l(\lambda_0 +h ) = l(\lambda_0)+h l'(\lambda_0) + \frac{1}{2}l''(\lambda_0)h^2 \]

Using the first and second derivatives above, we can easily make the second order Taylor polynomial for \(l\). taylor2 encodes this, where the argument \(x = \lambda_0 +h\).

Task

Complete the function below, which should return the value of the second order Taylor approximation to the log likelihood evaluated at \(\lambda_0\).

Overlay a plot of your approximant at \(\lambda_0 =1\), say, on a copy of the log likelihood curve given above.

taylor2 <- function(lambda0, x, y){

}

Solutions

taylor2 <- function(lambda0, x, y){
  n <- length(y)
  llk_pois(lambda0) + 
    (sum(y)/lambda0 -n)*(x-lambda0) -
    0.5*sum(y)*((x-lambda0)/lambda0)^2 
}

plot(lambda_seq,llk_curve, type= "l",xlab="lambda",ylab="log lk")
lines(lambda_seq,
      taylor2(1,lambda_seq,y), col = "red")

You should find that your approximant is close to the log likelihood for a small region around \(\lambda = 1\).

The idea behind Newton’s method is as follows. Even if we can’t maximize the log likelihood, we can always maximize the second order Taylor approximant - this is just a quadratic after all. We then take the maximum of our approximant as a new approximation to the maximum of the log likelihood, and repeat.

Task

Complete the function below, which returns the maximum point of the Taylor approximant - you should be able to determine easily where this occurs, in terms of the derivatives of the log likelihood. Add a dashed line showing the maximum on the plot above.

newton_pt <- function(lambda0, y){

}

Solutions

newton_pt <- function(lambda0, y){
  n<-length(y)
  lambda0 + (sum(y)/lambda0 - n)*lambda0^2/sum(y)
}

plot(lambda_seq, llk_curve,
     type = "l",
     xlab ="lambda",
     ylab ="log lk")
lines(lambda_seq,
      taylor2(1,lambda_seq,y),col="red")
abline(v = newton_pt(1,y),
       col="red",
       lty=2)

We might reasonably hope that this point, as the maximum of our approximant, is a better approximation to the maximum of the original function than the point we started with. Here, this is indeed the case.

Task

Do six (or so) steps of the process above: at each step, find the maximum of the Taylor approximant and use it to update your estimate of the maximum. Start the process at \(\lambda_0 = 0.4\), further away from the maximum, so that you can see what is going on; it converges quickly. You are aiming here to make a figure similar to Fig 2.4 of the notes.

Solutions

x0 <- 0.4

par(mfrow = c(2,3))
for(i in 1:6){
  plot(lambda_seq,
       llk_curve,
       type = "l",
       xlab = expression(lambda),
       ylab = "log lk")
  lines(lambda_seq,
        taylor2(x0, lambda_seq, y),
        col = i)
  x0 <- newton_pt(x0, y)
  abline(v = x0, col=i, lty=2)
}

3.3 Newton’s method for the Poisson GLM with log link

In a Poisson GLM with the log link, the expected response of the \(i\)th observation is

\[ E[Y_i] = \lambda_i = \exp(\eta_i) = \exp(X_i\beta). \]

As in the previous example, we can easily write down the log likelihood as a function of \(\beta = (\beta_1,\ldots,\beta_p)\),

\[l(\beta) = -\sum_{i=1}^n \lambda_i + \sum_{i=1}^n y_i \log(\lambda_i) -\sum_{i=1}^n \log(y_i!) .\] Newton’s method carries over nicely to the multivariate setting. In terms of the gradient and hessian of \(l\), the update step is just

\[ h = - (\nabla^2 l)^{-1} \nabla l \]

Task

For the Poisson GLM above, show that the gradient and hessian of \(l\) are given by

\[ \nabla l = X^T(y -\lambda), \qquad \nabla^2 l = - X^T W X, \]

for a suitable matrix \(W\).

Solutions

We can quickly implement Newton’s method with the following functions:

grad_pois<-function(beta,X,y){
  lambda <- exp(X%*%beta)
  t(X)%*%(y - lambda)
}

hess_pois<-function(beta,X,y){
    lambda <- as.vector(exp(X%*%beta))
    -t(X)%*%diag(lambda)%*%X
}

Here we simulate some data from the Poisson regression model,

n <- 50
x <- seq(0,1,length.out=n)
beta0 <- 0.5
beta1 <- 0.3
beta <- c(beta0, beta1)
X <- cbind(1,x)
eta <- X%*%beta
y <- rpois(n, lambda = exp(eta))

Task

Use the functions above to do ten iterations of Newton’s method starting at the point \((\beta_0, \beta_1) = (0, 1.5)\). Plot the trajectory on a contour plot of the log likelihood.

beta<-c(0, 1.5)
beta_store <- matrix(nrow=10,ncol=2)
for(i in 1:10){
  beta_store[i,]<- beta
  beta <- beta - solve(hess_pois(beta,X,y),
                      grad_pois(beta,X,y))
}

glm(y ~ x, family = "poisson") # R function for fitting GLM

## 
## Call:  glm(formula = y ~ x, family = "poisson")
## 
## Coefficients:
## (Intercept)            x  
##      0.4635       0.2434  
## 
## Degrees of Freedom: 49 Total (i.e. Null);  48 Residual
## Null Deviance:       57.94 
## Residual Deviance: 57.48     AIC: 169.4

llk_fun <- function(beta0,beta1) llk(beta0,beta1, x, y)
llk_vec <- Vectorize(llk_fun)
beta0_seq <- seq(-4,4,length.out=100)
beta1_seq <- seq(-2,2,length.out=100)
llk_contour <- outer(beta0_seq,
                     beta1_seq,
                     FUN = "llk_vec" )
contour(beta0_seq,
        beta1_seq,
        llk_contour,
        levels=seq(-130,-80,5),
        xlab = "beta0",
        ylab="beta1")
lines(beta_store[,1], beta_store[,2], col = "red")
points(beta_store[,1], beta_store[,2], col = "red")